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Abstract. Truncated Taylor series representations of invariant manifolds are 
abundant in numerical computations. We present an aposteriori method to 
compute the convergence radii and error estimates of analytic parametrisa- 
tions of non-resonant local invariant manifolds of a saddle of an analytic vector 
field, from such a truncated series. This enables us to obtain local enclosures, 
as well as existence results, for the invariant manifolds. 
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1. Introduction 

The invariant manifolds of a saddle of a vector field are very important objects 
for the understanding of the global dynamics of the flow generated by the vector 
field. The invariant manifolds divide the phase space into regions with different 
behaviour. The simplest picture is in the plane, where the invariant manifolds, 
the separatrices, of the saddles of the system, together with the limit cycles, de- 
compose the phase plane into connected components where the trajectories have 
similar a and ui limit sets. A standard reference on invariant manifold theory is 
In higher dimensions the structure of the invariant manifolds is, typically, 
much more complicated. The fundamental theorem about hyperbolic saddles is 
the stable (unstable) manifold theorem, see e.g. [Ill 120) . which states that locally 
at a hyperbolic fixed point there exist manifolds of dimensions d s and d u , denoting 
the number of negative and positive eigenvalues of the linearisation of the vector 
field at the fixed point, such that the tangent spaces of the stable and unstable 
manifolds at the fixed point are the negative and positive eigenspaces, respectively, 
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of the linearisation. In addition, these manifolds can locally be described as graphs 
of functions from the negative eigenspace to the positive eigenspace, and vice versa. 

To compute global invariant manifolds, one typically starts with an approx- 
imation of the local invariant manifolds lying in the corresponding eigenspace of 
the linearisation, and expand the global invariant manifolds step by step from the 
local one. For a review of a plethora of such methods see |13| . 

For obtaining approximations of (local) invariant manifolds there exists many 
references, e.g. [2J 01 03 [HI [S3J . Few methods exist, however, that can rigorously 
compute enclosures of the local invariant manifolds, which is our current objective. 
Some such methods are [TBI [19\ [26] . 

We present a method to compute the convergence radii together with explicit 
error estimates of the parametrisations of the invariant manifolds; for some meth- 
ods to compute such parametrisations see e.g. [4} [5] 123]. The parametrisations that 
we study are constructed such that the negative and positive eigenspaces of the 
linearisation at the fixed point are invariants of the flow of the vector field. This 
is a much weaker requirement than to completely linearise the vector field, as can 
be done, according to Siegel's theorem [3T], under certain Diophantine conditions 
on the eigenvalues. The idea to compute a close to identity transformation that 
removes all terms necessary for the transformed equation to have this property 
has appeared in [53]; in [35] the resulting vector field, after this close to iden- 
tity transformation, was named a robust normal form. These linearisations can 
be seen as a special case of the parametrisations of invariant manifolds in [4] [5J, 
where higher order conjugacies are also considered. The constructive method to 
compute convergence radii and error estimates presented in this paper, however, 
are much easier to implement and compute than the aposteriori convergence op- 
erator from [U[S]. Since the case of conjugacy with a linear flow on the invariant 
manifolds is probably the most common, the fast and simple results from this 
paper, that directly translate into an algorithm to compute convergence radii and 
error estimates, should potentially be very useful. 

This paper is organised as follows: in Section [2] we introduce the necessary 
notation, recall the necessary concepts about robust normal forms from [121 125] . 
and state our main result on the existence of analytic parametrisations. In Section 
[3J we prove the main theorem. The proof of the theorem is constructive and in 
Section [4] we describe an algorithm that implements the proof. Finally, in Section 
[5j we calculate the local invariant manifolds in a simple planar system, similar to 
the discrete system studied in |18 [ fl9 | [26]. 



2. Statement of the results 

Consider a vector field in R d of the following form: 

z = Az + F(z), (1) 

with A e S, where S :— {diag(Ad B , . . . , A 1? . . . , ^d u ) : Ad s < • • • < Ai < < 
Hi < ■ ■ ■ < Md u }i ancl where F is an analytic function, with F(z) = 0(z 2 ). Note 
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that any vector field with a saddle fixed point, with distinct real eigenvalues, can 
(locally) be brought into this form by an affine change of variables. We decompose 
z in the stable and unstable coordinates, z — (x, y) £ M. ds x M d " . We use = 
(0, . . . , 0, 1, 0, . . . , 0) to denote the vector in R d with its ith component equal to 1. 

The structure of the parametrisation of the invariant manifolds that we are 
computing is based on the close to identity change of parameters associated with 
the robust normal forms studied in |121 |2~41 125] . In order to simplify the formulae, 
we use vector and multi-index notation. In this section we revise and adapt the 
necessary notation and results from [12\ 124). and state our main result. 

The structure of ((TJ) implies that the stable and unstable manifolds at the 
origin are tangent to the coordinate axes. As discussed in the introduction, we seek 
a parametrisation of the stable and unstable manifolds, i.e., we want to compute 
maps 4> and ip: 

: R ds — ► R d 
ip-.W 1 " — > R d 

that are such that: 

W£ c = {(£,0) + <K£):£et/cR d *}, (2) 



W l l c = {(0, V ) + ^( V ):T 1 eVcR d "}. (3) 
Note, this means that the stable and unstable manifolds are not represented 
as graphs; we compute parametrisations of the invariant manifolds, i.e., we allow 
nonlinearities in the stable coordinates of <fi and the unstable coordinates of tp, 
respectively. 

We require that <f> — 0(f 2 ) and ip = 0(t] 2 ). The maps <fr and ip determine a 
close to identity change of coordinates in M. d : 

e(t,v) = ((,,v) + m + Hv)- (4) 

The idea of the parametrisation is that in (£, ^)-coordinates the local stable 
and unstable manifolds should be given by E s (= m. d °) and E u (= the stable 

and unstable tangent spaces at the fixed point. £ should be interpreted as the 
nominally stable, and r\ as the nominally unstable coordinates. Since (j) and ip do 
not have any constant or linear parts, the pullback of the original vector field using 
has the following form: 

8*(A + F) = A + G. (5) 
The formal power series for the nonlinear part of the vector field in the new 
coordinates is: 

oo 

G = 9mC m (6) 

|«i|=2 

In order for the local invariant manifolds to be of the forms ([2]) and ([3]) G must 
be of order 0(min(|£|, This means that if there exists i, 1 < i < d, such that 
g m Gi is a non-zero coefficient in the formal power series of G, then \m s \ > 1 and 
\m u \ > 1. 
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Thus, 



G\ 



and G\e 



We call the non-negative number \m\ 



(7) 

the order of to, and 



define the set N 2 = {to G N 2 



>2}. 



We split the space of multi-indices into the sets 



V := {m e N 2 
U := {me N 2 
where V is further decomposed into 
V s = {to G V : \m u \ = 0} 



E 



\m s \ = or |to„| = 0}, 
m s | > 1 and \m u \ > 1}, 
= V S UV„, where 
and V„ = {m G V : |m s | = 0}. 
We can now define the set of admissible linear parts of ([1} that we consider: 
T a := {A G S : to e V => to s A - A 4 ^ 0, 1 < i < 4}, 
J„ := {A 6 5 : m G V =>■ m„/i — /ii ^ 0, 1 < i < d u }, 
J 7 = T s r\J 7 u . 

We will often use the notion of filters of a (formal) power-series: if f(z) = 



m|>2 



a m z m , we define 



[/]c = 



meu 



[/]v = 



mev 



and [/] m = a r 



Note that in this notation we have 



^, and [G]o = G. 



Also, we let denote the partial sum of the terms of / up to order d. We 
use the norms \z\ = maxKKd {l^ill an< ^ ||/||r — max {|/( z )| : \ z \ < r }- The r-disc 
is denoted by *B r . If X is a set and r G R, we denote by rX the set {rx : x G X}. 
The smallest integer, n, larger than a real number, r, is denoted by n = |Y] . 

Let a — (ai, . . . , afc) G K. fc , we say that a is A- finitely rationally independent 
if aj ^ ma, for all multi-indices m such that 2 < |m| < A, and all 1 < i < k. 
Finally, let Q : Z + — >• R> be defined as 



fi(fc) := min(|fcAi - X ds \ , fc/ii - Hd u \) ■ 

We will use the following lemma, which essentially is a reformulation of | 
5.1]. 



(8) 
Lemma 



Lemma 2.1. Assume that A is 



Ai 



-finitely rationally independent and ji is 



finitely rationally independent. Then, A G T . Furthermore, for all multi-indices 
to G V with orders \m\ > max ^ 4^- , ^ we /iGwe the following sharp lower 

bound: 

for all v G {AJ U (9) 



|to • (X,fx) —v\> f2(|m|), 
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From this lemma it clearly follows that T is open. In addition T has full 
Lebesgue measure in 5, since it is constructed by removing countably many lines 
from S. 

We are now ready to state our main theorem : 



Theorem 2.2. Given a system z = Az + F(z), where F(z) = X)| m |>2 c " 



analytic function, A £ T , and a natural number n\ > max ( 
exists analytic parametrisations of the stable and unstable manifolds of the forms 
(0) and {3|) ; converging on the disk Q3 re , with 

[eiV 1+ 7, lei^ 1 



there 



2<\m 3 \<n 1 



-1 



iKfl)€ ]T Pm^+rerf) 1- M ) ><®i (11) 

2<|m„|<ni V r B/ V r © / 

/or a computable positive real number re ■ 

To prove the convergence of the change of variables we proceed as in e.g. 
[TUll22| . and use the method of majorants. If /, g : C d —> <C d are two formal power 
series, and \f m \ < g m for all multi-indices m, and all the coefficients of g are real 
and positive, we say that g majorises /, denoted by / -< g. Thus, the convergence 
radius of / is at least as large as that of g. We will majorise in two steps; given 
some / : C d — > C d , we construct g : C d — > C such that /j -< g, for all i, and then 
construct h : C — > C such that g(z, z) -t, h(z). 



3. Proof of the main theorem 

Let z = (x, y) be the original coordinates, and £ = (£, 77) the coordinates in the 
domain of 6. Recall, 0(£, rf) = (£, rj) + + il>(r}), where [4>]v s = <j>: an d [i>]v v — ^P- 
By inserting z = 0(C) into {!]), differentiating, and comparing the sides, we get: 

£0(C)C = i = A0(C)+F(0(C)). 

Inserting this expression into ([5]) yields: 

DQ(()A( + DQ(C)G(0 = A0(O + F(0(C)), 

we reorder the terms and get: 

Let La and K\ be the operators 

L A 0=[Ity(OAC-A#0]v., (13) 
K A i/j = [2tyfo)AC - A^(»?)]v„, (14) 

where we note that 

Lxir-ei) = (m s A- (A,/i) i )^ m «e i (15) 
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K A {r, m "ei) = (m uP - (A, ^) i )^ m "e l . (16) 

Thus, since is a series in £ ms terms, and ip is a series in r] mu terms, the left hand 
side of (TT2")) can be written as L\4> + K A ij). Furthermore, 

[L A + K A ip] v = 0. 
This means that we have to construct G such that 

[F(9(C)) - D0(C)G(C)]u = 0. (17) 

Recall, we want to compute a normal form ([3]) such that G = 0(min(|£|, \r)\)), i.e., 
[G]u = G. To be able to do this we have to construct <j> and ip such that 

L\4> + K A ^ = [F(9(C)) - D6(C)G(C)] V - [F(Q(C))]v (18) 

To be able to simplify (TT5)) by decoupling the various terms into groups, we note 
that V s and V u are invariant under F, i.e., if m G V s , then to = (rn s , 0), and for 
any i 

oo 

F{z m e i )=F{x m °e i ) = £ c„(z m » ei ) n , 

|n|=2 

thus 

[F(x m 'ei)] v , =F(x m "ei), 

and similarly for Y u . Therefore, by filtering on the component level, we get the 
following two functional equations for and ipf. 

(L A cj ) ) l = [F((tO) + <j ) (0)}v B (19) 

(K A ip)i = [Fi((0,r,)+ij(r,))} Yu . (20) 
It follows that G should solve 

G(C) - me(C))} V - D{4>{0 + ^(r?))G(C). (21) 

Since <f> and ip do not contain any linear terms, this means that the coefficients of 
the formal power series of G can be computed recursively as 

9m - [[F<(e(C))]c- D{cj>{t) + ^))G[H-i]( C) " . 

Note that and (|2"0)) are the same formulae as the linear case of formulae 
[4, Equations (3.5)-(3.8)], but we have included their derivation for completeness, 
and to set the notation. 

Since A G J, and [<p]y 3 = 4> and [ip]y u = ip by construction, we can solve 
(fH?)) and PD|) recursively. By computing 0, -0, and G using the recursive formulae 

(f2"0"|) , and (f2"Tj) . we get formal power series with the properties described in 
the introduction, i.e., 

e*(A + F)| Bi =diag(A dsJ ...,A 1 ,0,...,0)(£,0), (22) 

and 

6*(A + F)\e u = diag(0, . . . , 0, Wl . . . , ^ ) (0, 77) . (23) 
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Thus, the stable and unstable manifolds are given by the parmetrisations ([2]) 
and ([3]), respectively. 

To bound the solutions of (jn|J) and fl20Jl we proceed as in [TH [Ml [25] , and 
prove the convergence of the change of variables using majorants and induction. 
Let 

N = max 











) 


( 


Ai 


1 


Ml 





be the constant from Lemma HOI from which the explicit lower bound holds. Recall 
that we have assumed that n\ > N. 
Let 

OO OO 

<t>i(0 = a i,mA ma and ^fa)= J2 h,m u V mu 

|m a |=2 m u |=2 

be the sought change of variables, where the on,m s and /3i, mtt with \m s \, \m u \ < n\ 
can be computed with any method that solves (|19[) and (j!20)) . To majorise the 
functions <fi and 0> we construct two one-dimensional functions cj> and Put 

"fe = max{|a iim |} and (3 k = V" max {|A, m |}- 

|m|=fc |m|=fc 

We then define 



^dfc£j fc and j>(uj) = T"' /3fckA 



fc=2 fc=2 

The and -0 are majorants of and -0, respectively. Although the convergence 
of the parametrisations of the local stable and unstable manifolds can be proved 
separately, for simplicity of the exposition, we henceforth study their convergence 
simultaneously. Therefore, let 

7fc = Oik + Pk, 

and define the joint majorant 

oo 
k=2 

To calculate oti im and /3j, m , with \m\ = k, we use the operators La and K\ 
defined by and (fLT| . respectively. Their evaluation reduces by (|19p and (j2"0)) 
to the evaluation of /c-Taylor models of Fj ((£, 0) + 0(£)) and Fj ((0, ry) + ip{j])), 
respectively. The action of La and K\ on monomials are given by (|15p and (I16D . 
and yield the following formulae for on <m and /3j lOT , respectively: 

^((^o)+^- 1] (o)U ) 



Am s - (A,/i)j 

[^((0, ?? )+^ fc ' 1 ](x))] 
- (A,/i)< 



(25) 
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Ik < CM k , 



<k<m. (26) 



Note that the coefficients at a certain level only depend on the previous 
levels. The reason is that F does not contain constant or linear terms. This in 
turn allows for a recursive solution scheme of (TTTJ1) and (|2T))) . given by (pM)) and 
(J^SJ), respectively. 

If ?ii is sufficiently large, then the first n\ terms of <p and ip produce a good 
approximation of a majorant x, and we use this to determine an approximate 
radius of convergence for \. The validity of this radius of convergence will now be 
proved. As a first step we determine, using a least squares estimator, constants C 
and M such that 

ni 
. 2 . 

The reason why we only use the terms from I and onwards, is that I should 
be large enough to capture transient phenomena in the sizes of the coefficients of 
X, so that the estimate from ([21H) is a tight bound on the coefficients in the tail 
of the power series of x- The least squares estimation is done in two steps: first 
a standard least squares approximation is computed, then we assume that M has 
been well approximated and increase C until (f2"6")l holds. Thus, a candidate radius 
of convergence is 

re = ± (27) 

which needs to be verified. 

We will consider a slightly larger majorant of %. If 

oo 

F( Z ) = CmZ" 1 , 

\rn\=2 



we define 
and set 



Emax {|ci, m |} , 
l<i<d 

\m\=k 



F{lu) = ^c fe uA 

fc=2 

F is clearly a majorant of Fi(z, . . . , z). Recall the definition of (|5J|, and let 

n = min min J l -, . 28) 

\2<|m|<Ar,i/e{A,}u{M,}. \m\ N J 

Note that is monotonically increasing for k > N, and |m • (A,/i) — v\ > fl(k) 
for all \m\ > N and ^ £ {A^} U {/ii}. Hence, 

\m ■ (A,/i) — > f2|?n|, 

for all |m| > 2, and all ^ € {A^} U {/J,i}- We recursively define a majorant 

oo 

c(cj) = 

fc=2 
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of xM by: 



Sk= <k 



Jk-i] 



M) 



for k > 2. 



In our proof of convergence of the parametrisations we will use a quadratic 
bound on F. If the convergence radius of F is s, we choose two other radii < 
s" < s' < s, and use Cauchy-type estimates on the p-tail of F on 23 s < valid on 
Q5 S ", and then require that 2re < s". 

Indeed, let iV<j(fc) denote the number of d-dimensional multiindices with ab- 
solute value k, then for any i and m, 



< 



l/ll 



(s')l m l 
and hence 

Therefore, we can bound the p-tail of F on s" as follows, 



< (^Er=, +1 w(^) 

We denote the bound on the tail 



fc-2 



v ' k=p+l v 



fc-2 



(29) 



and define: 



fc=2 

Clearly, \F(z)\ < A s »\z\ 2 , on <B S ». 

If possible put s" = 2re, otherwise take s" as large as possible and put 
re = s" 12. If re > xx^ — , decrease re until 

' e £ nb- <31) 

To prove the convergence of a we proceed as in jTOJ [13 ■ The definition of the 
5k 's imply that (formally) the following equation holds: 

oo oo 

Slk8 k uj k = Y Cfe + a(oj)) k , 

k=2 k=2 

where we note that the left hand side is equal to fluja' (ui). Hence, a satisfies the 
following differential equation 
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, F(u) + ct(lj)) , , 

= no, CT(0) = °- 

The fact that neither a nor i* 1 have any constant or linear parts imply that 
the following inequalities holds for any partial sum (the first is an inequality since 
the right hand side, in general, also includes some higher order terms): 



,</i[.)<^""'"" 



and 

0<aW(u)<uo'W{w), 0<uj. 
Together they imply that 

We will now use our quadratic bound on F. Assume that cr[ fe -i]( re ) < re, 
for some k (this trivially holds for k = 2), then 



F(r e + cT^jre)) 

n 

F(r e + r e ) 
Q 

A 2rB (r e +r ) 2 . 

n - re ' 

where the last inequality is due to (|31j) . Hence, by induction, cr(re) < re, and since 
all the coefficients 8k are positive this implies that a is analytic on *8 re . It follows 
from the convergence of cr, by tracing the sequence of majorisations backwards, 
that is analytic on 58 re . 

The remainder terms in Theorem 12.21 are found by using Cauchy bounds on 
cr. Since <j(lj) < r& on Q5 re , the convergence radius of a is larger than re, and we 
have that 

Since we use the supremum norm, the uncertainties can appear in any component 
simultaneously. Therefore, the remainders are added as 25 1 scaled with geometric 
bounds given by the bound on the growth of the 5^s. 



< 

< 

< 
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4. Algorithmic aspects 

The main application of the convergence proof given in this paper is that it is con- 
structive and suitable for implementation on a digital computer. We summarise 
the key points of the proof given in the last section and compile it into an algo- 
rithm, computing re from the formulation of Theorem 12.21 Since re tends to be 
slightly pessimistic, the algorithm also computes the constants C and M. They 
are candidates for the geometric bound on the tail of O. In general re < -h- This 
algorithm is given as Algorithm [1] 

5. Example 

An algorithm proving the conditions of Lemma |2~T1 and Theorem 1 2 . 21 has been im- 
plemented in a C++ program using the C-XSC package [5J[9] for interval arithmetic 
[TJ [T5l [16l [17] . For automatic differentiation [8] we use a modified version of the 
Taylor arithmetic package [3]. 

There are several methods to compute local invariant manifolds of discrete 
dynamical systems, see e.g. [TH1 HH US]- In principle, these methods can also be 
used for continuous dynamical systems by studying the time-i map of the flow for 
some t. To compare with the method in [26J, which is also able to treat flows, we 
study a vector field of the same form as the discrete dynamical system studied in 



Using ri\ = 81, we compute (the computation takes a few seconds) re = 
0.023, M = 2.69 and C = 0.30. These values yield the following bound on the 
error terms in Theorem 12.21 



The image OCB0.023) is shown in Figure [T] The image of 0(*Bo.37) given by 
the candidate radius of convergence jr is given in Figure [21 By inspection we see 
that the image 0(*Bo.o23) contains the ball So. 02; this can be compared with the 
convergence radius 0.18 with the bound 0.241138 on the Lipschitz constant [2"7] 
for the cone enclosures of the local invariant manifolds using the method from 
|26j . This indicates that a method to enclose local invariant manifolds on a larger 
domain could be to use the method [26] outside of the 9-image of the result of our 
method. 
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Algorithm 1: Implementation of the proof of the main Theorem 
Data: A, F, n\ 

Result: <t> [ni \ V [ni1 , re, C, M 
l for k — 2 to m do 
2 
3 



for i — 1 to d do 

forall |m a = k do 

[F 1 ((5,o)+^[ fc - 1 l(0)] r 



end 

forall |m„| = k do 

[F,((0,»,)+V Ifc - 11 (^))] r 



end 



end 



7fc = E|m s | = fc maX l<*<4K™J} +E|m u |=fc m aXl<i<d{|/3i,m u |}. 



r 

8 
9 
10 

11 end 

12 n = |m/2j; 

13 for k — no + 1 to m do 

14 | B{k,l) = 1, S(fe,2) = k, b(k) = log 7fe 

15 end 

16 (log C, log M) = (B T B)- 1 b ■ 

17 for fc = no + 1 to ni do 

18 if 7 fc > C*M fc then 

19 | C* = 7fc/M fc 

20 end 

21 end 

22 for k — 2 to p do 



Ck 



■i | m| — fc 



maxi<i< d {]ci, m |} 



24 end 



25 



fi = mm (min 2 <| m | <JViI , e{Ai}u{w} . Iffi , 



26 r e = ± 

27 repeat 



28 

29 

30 
31 
32 
33 
34 
35 



Compute A p 2rei ; 



1 Converges = True 



AS. 



n then 



else 



Converges — False; 
re = 0.95re; 



end 



36 until Converges = True ; 
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Figure 1 . Enclosures of the local stable (blue) and unstable (red) 
manifolds in the example. As the local invariant manifolds ap- 
proach the converges radius of O, the uncertainty of its location, 
given by the remainder term in Theorem l2.21 becomes unbounded. 
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Figure 2. The local stable (blue), and unstable (red) manifolds 
on their heuristic domain of existence. 
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